      subroutine iexpso(iin, iout, popp, popq, ioppq, iexppq,  
     &   popr, pops, ioprs, iexprs, sym) 
C     this program is to do the expand of I(pqrs) with iop 
C     if iexppq = 1, then expand p<q or p<=q to p,q  
C     if iexprs = 1, then expand r<s or r<=s to r,s  
C-----attention! this program can only deal with on side. 
C-----for example: iexppq/=0, then iexprs must equal to 0 
C-------------------
      use mod_symm
      use mod_iop
      implicit none 
      INTEGER, external :: IRPDSO
C----------------------------------

      integer, intent(in) :: popp(8), popq(8), popr(8), pops(8) 
      integer, intent(in) :: ioppq, iexppq, ioprs, iexprs 
      integer, intent(in) :: sym
      real*4, intent(in),managed :: iin(*)
      real*4, intent(out),managed :: iout(*)
      DOUBLE PRECISION, EXTERNAL :: DNRM2
C
      integer irrepr, irrepl, i, numl, numr, numl2, numr2, popl
      integer nsize1, nsize2, ioff, ioff1, nsize
C
      if(ioppq/=0.and.ioppq/=1.and.ioppq/=2) then
         write(6,*)'error of ioppq in isymtrso'
         stop
      endif
      if(ioprs/=0.and.ioprs/=1.and.ioprs/=2) then
         write(6,*)'error of ioprs in isymtrso'
         stop
      endif
C
      if(iexppq==1.and.ioppq==0.and.iexprs==0) then 
      ioff = 1 
      ioff1 = 1
      do irrepr = 1, nirrep 
         irrepl = dirprd(irrepr, sym) 
         numl = irpdso(irrepl, popp, popq, ioppq) 
         numr = irpdso(irrepr, popr, pops, ioprs)
         numl2 = irpdso(irrepl, popp, popq, iexppq) 
         call symexp2so(irrepl, popp, numl2, numl, numr, 
     &              iout(ioff1),iin(ioff))  
         ioff = ioff + numl*numr*isd 
         ioff1 = ioff1 + numl2*numr*isd
      enddo
C
      elseif(ioprs==0.and.iexprs==1.and.iexppq==0) then
      ioff = 1 
      ioff1 = 1
      do irrepr = 1, nirrep 
         irrepl = dirprd(irrepr, sym) 
         numl = irpdso(irrepl, popp, popq, ioppq) 
         numr = irpdso(irrepr, popr, pops, ioprs)
         numr2 = irpdso(irrepr, popr, pops, iexprs) 
         call symexpso(irrepr, popr, numl, iout(ioff1), iin(ioff))
         ioff = ioff + numl*numr*isd
         ioff1 = ioff1 + numl*numr2*isd
      enddo
C
      elseif(iexppq==1.and.ioppq==2.and.iexprs==0) then
      ioff = 1
      ioff1 = 1
      do irrepr = 1, nirrep
         irrepl = dirprd(irrepr, sym)
         numr = irpdso(irrepr, popr, pops, ioprs)
         numl = irpdso(irrepl, popp, popq, ioppq)
         numl2 = irpdso(irrepl, popp, popq, iexppq) 
         call exp4symso(irrepl, popp, numl, numl2, numr,
     &              iin(ioff),iout(ioff1))
         ioff = ioff + numl*numr*isd 
         ioff1 = ioff1 + numl2*numr*isd
      enddo
C
      elseif(ioprs==2.and.iexprs==1.and.iexppq==0) then
      ioff = 1
      ioff1 = 1
      do irrepr = 1, nirrep
         irrepl = dirprd(irrepr, sym)
         numl = irpdso(irrepl, popp, popq, ioppq)
         numr = irpdso(irrepr, popr, pops, ioprs)
         numr2 = irpdso(irrepr, popr, pops, iexprs) 
         call exp3symso(irrepr, popr, numl, iin(ioff), iout(ioff1))
         ioff = ioff + numl*numr*isd
         ioff1 = ioff1 + numl*numr2*isd
      enddo
C
      else
         write(6,*)'attention!!!!!, error in iexpso'
         stop
      endif
      return 
      end  
